Multistability in Networks of 
Weakly Coupled Bistable Units 



R.S. MacKay^ and J. -A. Sepulchre 



a,b 



^ Nonlinear Systems Laboratory, Mathematics Institute, University of Warwick, 
VO ' CV4 7AL Coventry, United Kingdom 

I ^ previous address: Universite Libre de Bruxelles, CP 231, Boulevard du 

Triomphe, B-1050 Bruxelles, Belgium 

c ■ 

O ■ Abstract 
(N , 

We study the stationary states of networks consisting of weakly coupled bistable 
^ , units. We prove the existence of a high multiplicity of stable steady states in net- 

\ works with very general inter-unit dynamics. We present a method for estimating 

• the critical coupling strength below which these stationary states persist in the net- 

work. In some cases, the presence of time-independent localized states in the system 
\ can be regarded as a 'propagation failure' phenomenon. We analyse this type of be- 

' haviour in the case of diffusive networks whose elements are described by one or 

two variables and give concrete examples. 

O 
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1 Introduction 



One of the simplest behaviours in nature which reflects underlying nonlinear 
mechanisms, is the phenomenon of bistability, that is the possibility for a 
system to present simultaneously two stationary states which are stable for 
the same set of system parameters. This type of behaviour has been identified 
in a great number of natural systems, ranging from laser physics to cellular 
biology. 

In physics, examples of bistable phenomena are the behaviour of a metal- 
lic strip under longitudinal squeeze, the input-output response of diode cir- 
cuits [1], the onset of electrical filaments in gas-discharge experiments [2], and 
the phase separation dynamics in condensed matter [3]. Another interesting 
example of bistability is studied in nonlinear optics [4] . For example, if a co- 
herent electrical field is injected into a ring cavity containing an absorbing 
medium, the latter can respond in two different ways for the same amplitude 
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of the driving field. In tlie first regime, tlie energy wliich is transmitted into 
the cavity is nearly totally reflected, whereas the second regime is character- 
ized by a high transmission rate. Such devices have been proposed for building 
optical memories. 

Other phenomena of bistability have been demonstrated in chemical systems. 
For instance, the bistability of the chloride-iodite reaction has been studied 
over a wide range of external constraints [5]. In biology also, interesting cases of 
bistability have been identified [6], with, for example, the membrane potential 
of neural cells which can sustain, under some circumstances, two stable equi- 
libria. This latter situation can be described simply by the FitzHugh-Nagumo 
equations [6]. 

In most of these systems, the mathematical description of the bistability phe- 
nomena is fairly well understood, and it involves typically a nonlinearity of the 
'cubic-type'. What is still under investigation, however, is the understanding of 
the possible behaviour that bistable elements exhibit, when coupled together 
in a large network. In this paper, we concentrate on the behaviour which ap- 
pears in a network of coupled bistable units, for small value of the coupling 
strength between units. Using the implicit function theorem we show that, 
for a wide class of interaction between network units, the individual bistable 
behaviour of the units persists in the system, if these units are weakly coupled. 
It gives rise to a high multistahility of the network stationary states. Some of 
these are very localized, for example, the continuations of the states for which 
all the units are in the same equilibrium, except for one unit which is in the 
other equilibrium. 

Many of the steady states of the network may disappear as the coupling 
strength is increased. Numerical experiments typically show that in this case, 
one can observe the onset of an activity front which propagates in the net- 
work. On the other hand, when the coupling strength is below a critical value, 
propagation becomes impossible, and we suggest that the presence of steady 
localized structures can be interpreted as propagation failure. This situation 
has been studied recently by several authors in the case where the interac- 
tion between units of the network is described by discrete diffusion [8]-[12]. 
For example. Keener [8] studied the propagation failure which can occur in 
a diffusive chain of excitable cells. He proposed a method for obtaining an 
analytical estimate of the critical coupling strength below which propagation 
is blocked. This method is valid for one- dimensional translation invariant dis- 
crete reaction-diffusion systems with one variable per unit. In the present 
paper, we propose a technique which allows one to extend this kind of estima- 
tion to networks with more general connectivity, not necessarily translation 
invariant and whose units are described by more than one variable. 

The paper is organized as follows. In section 2, we formulate the problem 
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as one of continuation of steady states from the uncoupled case. Section 3 
proposes a general scheme for estimating a critical coupling strength up to 
which all the steady states persist. Section 5 applies this scheme to networks 
described by two variables per unit interacting by means of discrete diffusion. 
In section 6, the same technique is applied for studying the persistence of 
steady states in networks consisting of bistable optical devices. In section 7 
we discuss the stability-type of the resulting steady states. Finally, section 8 
is devoted to conclusions. 



2 Continuation from the uncoupled limit 



In this section we show how steady states of an uncoupled network of bistable 
units can be continued to weak coupling, for a range of coupling which depends 
only on bounds on the dynamics of the individual units and the norm of the 
coupling, regarded as an operator on a sup-norm Banach space. We emphasize 
that the size of the network does not enter explicitly. 

We consider a network consisting of a large and possibly infinite number of 
units coupled together. The units of the network are indexed by an index set 
S (e.g. the integers Z or Z^), and the state of the unit n e S is described by a 
vector Xn = {vn, • • • , Wn) G R"^. Then, a state of the whole network is defined 
by an indexed set X = {xn)neS of m-dimensional vectors. We choose a norm 
I • I on R"*, and extend it to the set of bounded states of the network by 

\X\ — sup \Xn\ . 
n 



With this norm, the set of bounded states of the network forms a Banach 
space that we denote by R"*^. Now we define the dynamics of the network. 

First we assume that the dynamics of each unit taken separately is described 
by the differential equation 

^ = fM , (2.1) 



where the function / : R™ R™ has the property of bistability. This 

means that the equation f{xn) = has two distinct solutions 

0, 1) such that all eigenvalues of the matrix D/(a;*^*^) have negative real part. 

Next we define the map F : R*"^ R"*® by 

F{X) = ifixn)Us . 
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Then we assume that the dynamics of the network of bistable units obeys a 
differential equation of the type 

d X 

-^ = F{X) + aK{X), (2.2) 



where a is a real parameter called the coupling strength and K : R™ — >• 
R™^ is a function which describes the interaction between units. This 
interaction function is possibly nonlinear. An example for the function K, of 
special interest in physics, is the discrete diffusion operator in a d-dimensional 
cubic lattice. In that case, each component of X is a linear function 

2d 

K„{X)^Dj2{xn^-Xn) , (2.3) 

k=l 



where D is a, m x m diagonal matrix specifying the 'diffusion coefficients', 
and the subscripts k in eq. (2.3) run over the 2d nearest neighbours of the 
unit ,T„ in the lattice. For d — 1, it gives the well-known finite difference 

Our analysis starts from the uncoupled case, a = 0. This provides explicitly a 
large number of stationary states. Indeed, for any choice of binary indexed set 
(4)fceS, ^fe = or 1, there exists a corresponding stationary state of eq. (2.2) 
given by 

Xo = {x^'^\es ■ (2.4) 



It satisfies the equation F{Xq) — 0. 

Our idea is that all of these steady states can be continued to steady states 
for small a. The strategy is borrowed from [7]. 

Theorem 2.1 Given a system of the form (2.2), 3 > such that any 
stationary state Xq of the uncoupled case has a locally unique continuation for 
all \a\ < ckq. 



Proof. This theorem is proved by a simple application of the implicit function 
theorem. Let us rewrite the right-hand side of eq. (2.2) as a map G : R"*^ x 
R R"^S defined by 

G{X,a) = F{X) + aK{X) . (2.5) 



The stationary state Xq is such that G{Xq, 0) = 0. The zero {Xq, 0) of G has a 
locally unique continuation around o; = if, by the implicit function theorem. 
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the derivative of G with respect to X, denoted by DG, is invertible at {Xq, 0) 
(with bounded inverse) and the map G is continuously differentiable. 

Now, DG'(Xo, 0) = DF(Xo) is block diagonal with blocks \ = D/(a;(^'=)), and 
each block is invertible since all the eigenvalues of Ji^ have strictly negative 
real part. Hence DG^^(Xo,0) exists and is block-diagonal, with blocks equal 
to J~^. Since the norms of the blocks are bounded, DG^^ is bounded. F is G^, 
since f is G^ . K is G^ by assumption. Hence G is G^. Therefore there exists 
an interval of values of a around such that the equation G{X, a) = defines 
a function X{a). 

A common interval [— ao, olq] can be found which works for all choices of steady 
state Xq because the estimates in the imphcit function theorem can be chosen 
independently of the choice ol Xq. □ 

Remcirk 2.2 The function X{a) is differentiable with respect to a and 

= -DG-VX, a) KiX) , (2.6) 

da 

as long as DG(X, a) is invertible. In the next section we will estimate a lower 
bound on a for which DG~^ exists. 

Remcirk 2.3 The result of theorem 2.1 is independent of the specific form of 
the interaction function K{X) (provided it is G^). Consequently, the theorem 

is valid for a class of network interactions which is much wider than the discrete 
diffusion defined by eq. (2.3). Note that when K is linear, the G^ condition is 
just that the sum of the norms of the elements along the rows of the matrix 
K be bounded. 

Remark 2.4 Theorem 2.1 remains true even if the local dynamics / varies 
a little from unit to unit. All that is required is uniform bounds on D/ and 
(D/)~^ in a neighbourhood of the equilibria of /. 



3 Procedure for estimation of ckq 

Theorem 2.1 means that in a network of weakly coupled bistable units, there 
exist a great number of nonuniform states of the system which do not evolve 
in time. Stability is proved in section 7. 

The persistence of all the stationary states will occur as long as the operator 
DG(X, a) remains invertible. The goal of this section is to estimate a value ao 
of the coupling parameter up to which all the stable steady states persist. It is 
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possible, and indeed likely, that they will persist further than oq. Furthermore, 
the true threshold in a for continuation of a steady state will in general vary 
from steady state to steady state. So should be interpreted just as a lower 
bound on the set of critical a's for all the steady states. 

The method is based on the estimation of the norms of each sides of eq. (2.6). 
Let us define the variable 

C^\X-Xo\, (3.1) 



If II DC ^{X, a) II can be bounded by a function P{(, a) and || K{X) 
function (5(C), then the solution Xq can be continued with 



by a 



da 



<p(C,a)g(c), 



as long as P{C, a) < oo. Let cto be the first a such that P becomes infinite 
along the solution of 

^ = ^(C,«)Q(C), (3.2) 



starting from ^(0) = 0. Then the solution X{a) can be continued at least up 

to CCq. 

To implement this procedure, we first note a useful bound for || DG~^ \\. By 
definition of the function G we have 

DG^DF + aDK . (3.3) 



So if DF is invertible and jaj || DK \\<\\ DF ^ || ^ then DG is invertible and 

\\^G~^\\<u , „ , , , „ 77. (3.4) 

" DF-i -1 - a \\DK\\ ^ ' 



This is easy to use because DF is block diagonal, so || DF~^ ||~^ is the minimum 
of the quantities || T)f{xn)~^ ||~^. Sometimes it is useful to obtain a shghtly 
stronger bound on T)G~^ by incorporating the diagonal part of Y)K into DF, 
i.e., write 

DX = Lo + Li , (3.5) 



where Lq is the block diagonal (i.e. on-site) part of the operator Y)K and 
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Li = DK — Lq is off-diagonal (i.e. inter-site part). Tfien tfie operator DG is 
decomposed as 

DC = A(a) + a Li , (3.6) 



wliere A (a) = DF + aLo is block-diagonal by construction. Thus the inverse 
A~^ is also block-diagonal, with blocks equal to (D/(a;„) -|- Q;(Lo)n)~^ that we 
assume to be bounded because this is true for small a by our hypothesis on 
the spectrum of D/(xW). Then, provided \a\ \\ Li \\<\\ ||~^, the norm of 
DG~^ can be bounded by 

||DG'-^||< „ I I II r II • (3-7) 

" " A-i(a) -1 - a Li ^ ^ 



After this, there are many ways of proceeding the bounds || A~^(a) ||~^, ||Li || 
and \K\. If the system has special features, one can obtain nice bounds on 
these quantities and then solve eq. (3.2) explicitly. We will give an example in 
section 4. If not, one can nevertheless always obtain bounds of the following 
form, for ^ < some Ci, since F and K are C^: 



A ^{oi)\\ ^ > ao — ai(^ — bi\a\ 
II Li{X) \\<bo, 
\K{X)\<co + c,C, 



(3.8) 
(3.9) 
(3.10) 



with gq, ai, bo, bi, Cq and Ci real numbers, Gq, bo, Cq > 0. Therefore, considering 
the norm of each side of eq. (2.6), and combining eq. (3.7)-(3.10), we obtain 
the differential inequality 



da 



< 



ao — ai( — b\a\ 



(3.11) 



with b — bo + bi. Starting with the initial condition ^(0) = 0, eq. (3.11) holds 
as long as the denominator remains positive and ( < (i. The case of equality 
in eq. (3.11) can be integrated by introducing a new variable s in order to 
separate the numerator and denominator of the fraction (i.e., replace by 
-^/-^). The integration of the resulting equation permits one to determine 
the point (^o, ckq) for which the denominator of eq. (3.11) vanishes for the first 
time. If we assume Co < Ci) we obtain 

b bci \ aiCo J 
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It gives a lower bound on a up to which DG(X(a),a) remains invertible. 
Using the binomial expansion, it is instructive to expand cvq as the series: 




(26 + ci) al 
6 afco 




) 



(3.12) 



2 oiCo 



We note that the parameters h and ci enter only from the second term. 

To end this section we notice that the solution of eq. (2.6) is not affected by the 
multiplication of DG and i^' by a common regular matrix M . Nonetheless the 
relation (3.11), which is in fact derived from eq. (2.6), can be altered by such 
a matrix multiplication. This property can be used in order to optimize the 
estimation represented by the inequality (3.11). In the next section we apply 
the scheme presented above in the case of bistable elements described by one 
variable and when the interaction between units is described by a diffusion-like 
operator, as defined by eq. (2.3). 



4 One variable systems 

To demonstrate the method, we begin by treating the simplest case: one- 
dimensional translation-invariant nearest-neighbour-coupled one-variable sys- 
tems: 

Vn = f{Vn) + Oi {Vn+l + Vn-1 " '2Vn) , (4.1) 

with i>„ e R and n E Z. We restrict attention to a > as this is the case 
of most physical interest. We suppose / has (at least) two stable equilibria, 
which can be taken without lost of generality at v^'^^ = and = 1 (with 
/'(0),/'(l) < 0). Let Iq.Ii be the largest intervals about 0, respectively 1, on 
which f'{v) < 0. We illustrate (3.7), and the special form of this problem will 
allow us to find good estimates for (3.8)-(3.10). 

The blocks of A are f'{vn) — 2q;. So, if Vn £ U /i, for all n, and if a > then 



inf|/'(^;,)|+2a 



Now II Li II = 2. Thus eq. (3.7) gives 



DG-\v,a)\\< 



1 



(4.2) 



mf.J/'(t'n)| ' 



for Vn e IqUIi , 
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and for a > 0. Note the cancellation of the "2a" s which is convenient, though 
not essential. Next, in the present case Kn{v) = (f„+i + Vn-i — 2v„) and 

\K{v)\<2 + iC, 



with ( defined by (3.1). We can do slightly better than this, however, because 
all equilibrium states for a > satisfy 

f{v) < , f{v) > , (4.3) 



where v = inf w„ and v — supVn- To prove this, note that the equation for an 
equilibrium state is 

fiVn) = -a {Vn+1 + Vn-1 - 2Vn) ■ 



Now, the relations (4.3) and the negative slope of / at f = and 1 imply that 
the equilibrium states obtained by continuation from v„ e {0, 1} (for all n) 
satisfy 

Vn e [0, 1] for all n . 



Hence along continuation we have 

\K{v)\<2. (4.4) 



Combining eqs. (4.2)-(4.4) we deduce that the continuation can be performed 
with 



with the function D{Q defined as 
D{C)^ inf \f'{v)\, 



as long as the denominator D{Q remains positive. The inequality (4.5) is easy 
to solve and thus continuation is possible at least as long as a < ao with 



Co 

i 



1 

2 / ^(0 dC , (4.6) 



2 





9 



where (^o is the first value of ( sTich that D{() = 0. Equivalently, is the 
distance from {0, 1} to the nearest critical point c of /. 



This result holds whether f'{v) is monotonous or not. If f'{v) is monotonous 
on [0, Co] and on [1-Co, 1], then D{C) = min(|/'(C)|, | /'(I -()!)• Different cases 
may be distinguished. The simplest case is when c < | and |/'(C)I ^ l/'(l~C)l 
for all C e [0,c]. Then 



Co c 

Jd{C) dc = / -/'(C) dC = -/(c) . 



A similar case with c > 1/2 is treated analogously. Thus, in this simplest 
situation, continuation is possible for 

a<ao= { ^ 2 4.7 

ifoi 



In other words, our estimate of ckq is the half height of the nearest critical 
point of / to or 1. 

Another example is to suppose c < | but 

|/'(C)l>l/'(i-C)|forCe [0,Cm] , (4.8) 
l/'(C)l<l/'(i-C)|forCe[C„,c]. 

Then, integral of -D(C) can be calculated again and (4.6) gives 

«o = ^ (-/(c) + /(CJ + /(I - Cm)) ■ (4.9) 



We will see an example of that type in section 6. 

Example 4.1 A typical example of a simple model of one- variable bistable 
dynamics is the cubic Nagumo model [6]. In this case / can be defined as a 
cubic polynomial p{v) of the form 

p{v) ^ v{v - r)){l - v), with < 77 < 1 . (4-10) 



The equilibria and 1 of p are stable as p'{0) = —f] < and p'{l) = 1 — rj < 
0. We restrict attention to < 77 < | (the complementary case is treated 
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similarly). Then the nearest critical point of p to {0, 1} is 



civ) = ^ (r? + 1 - ^l + ?7(?7-l)) - I ■ 

Furthermore it can be shown that in the present case |p'(C)| < — C)l fo^' 
C < |. Hence, the result expressed by eq. (4.7) applies and the continuation 
from the uncoupled limit of cubic Nagumo units is possible at least for a < 
ao{r]) = -HM^. The analytical expression of q;o(?7) is rather complicated but 
we can expand it as a power series in rj: 

«o = ^-?^ + C?(77^). (4.11) 



This value of ao is in agreement, at the lowest order, with the value rj^/S 
predicted by Keener with a different approach [8]. The method proposed by 
Keener relies on a Smale's horseshoe type construction. Although result (4.11) 
is not quite as good as Keener's [8], our method generalises easily to many sys- 
tems inacessible to Keener's method. For example, take S to be a translation 
invariant lattice where each site has N neighbours with discrete diffusion 

N 
3=1 

where rij are the neighbouring sites to n (e.g. N — 2d for a d-dimensional cubic 
lattice) . Then the only change required to the above analysis is to replace the 
"2" in eqs. (4.4)- (4.6) by "N". 

Extensions to one-variable systems on non-translation invariant lattices or 
with other forms of coupling are also relatively straightforward. For example, 
take S to be any graph, and 

3 

with {Dnj > 0; n,j G S} a set of real numbers characterising the coupling 
between the bistable units of the system. As a generalisation of eq. (4.1), we 
consider the network of coupled bistable units described by the system: 

Vn = f{Vn) + aKniv) , 
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with a > 0. Then, doing the same reasoning as for eqs. (4.2)-(4.5), we obtain 
the following differential inequality: 



da 



< 



M 



D{C) - (M - m)a ' 



(4.12) 



where M — max„ Dnj and m — min„ J2j Dnj. The case of equality is easily 
integrated and we can express ccq under the form 



where is the first value of ^ for which the denominator of (4.12) vanishes. 
Extensions to networks where the local dynamics Vn = fn{vn) are not exactly 
the same at each site are also feasible. 



5 Diffusive networks with two-variable units 



The aim of this section is to apply the differential inequality of the type (3.11) 
to networks of bistable elements described by two variables per unit. Thus we 
suppose that the vector Xn = {vn, Wn) G and the function /(x„) introduced 
by eq. (2.1) has two components f{xn) = {h{xn), g{xn)). As the two stationary 
states x^'^\ {i = 0, 1) of eq. (2.1) are supposed to be stable, the product of the 
two eigenvalues of D/(x'^*)) is positive, and thus: 

>1« =detD/(x®) > 0. (5.1) 



On the other hand, wc suppose here that the units of the network are coupled 
with the discrete diffusion operator defined by eq. (2.3) in a rf-dimensional 
lattice, and we assume that the diffusion coefficients of the variables Vn and 
Wn are in the ratio 1: 6, respectively. 

In order to fix the ideas, we rewrite eq. (2.2) with the present hypotheses, and 
in the case of o? = 1, i.e., for a chain of bistable units: 



dv 

= h{Vn, Wn) + a {Vn+1 - 2Vn + Vn^l) , 

= g{vn: Wn) + aS (w„+i - 2w„ + W„_i) , (5.2) 
with n e Z . 
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Our analysis is valid, however, for arbitrary d. 

Example 5.1 We consider the FitzHugh-Nagumo model [6] which is a simple 
two- variable model of bistable dynamics which generalises example 4.1. In this 
case the functions h and g can be defined as 



where p{v) is a 'cubic-shape' function, e.g., p{0) = pij]) = p{l) = 0, p{v) < 
for < V < 7], and p{v) > for f] < v < 1. The parameter e characterizes the 
time-scale difference in the dynamics of variables Vn and Wn, and the parameter 
7 is chosen such that the curves defined by h{v, w) — and g{v, w) — in the 
plane {v, w) have three distinct intersections. 

The eq. (5.2)-(5.3), with 6 = and e < 1, are used in theoretical biology 
for modelling pulse propagation along the membranes of nerve cells [6]. Note 
that when 6 = 0, the study of the equilibrium states of (5.2) reduces to a 
one-variable problem, as Wn can be eliminated using g{vn, Wn) = 0. So to 
have genuinely two- variable problem we take 6^0. This is relevant to several 
interesting problems. 

Now let us take a binary indexed set {ik)ki^s-, {ik = or 1), and the corre- 
sponding stationary state Xq [eq. (2.4)] of the network defined by eq. (5.2) 
in the uncoupled hmit a = 0. By theorem 2.1, there exists a locally unique 
continuation X{a) of Xq for a small. We will estimate a lower bound on the 
coupling parameter a for which the stationary state X{a) persists. We will 
restrict attention to a > as this is the case of main physical interest. 

Following the procedure described in section 3, we denote the right-hand side 
of eq. (5.2) by the function G{X,a). Then we decompose DC according to 
eq. (3.3) and we denote DF by A. Here, each block of the operator A is a 2 x 2 
matrix defined as 



In this expression, the subscripts v and w denote partial derivatives, e.g., 
hy = |^(x„). If each block A„ is invertible, so is A and we can write 



h{v, w) = p{v) — w , g{v, w) = e{v — ■jw) , 



(5.3) 




(5.4) 




(5.5) 
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where 5*0 and form the partition of S induced by the indexed set {ik}keS, 
i.e., /c e 5*0 if ifc = and k & Si if ik — 1. 

We are free to choose the norm in the two-dimensional state space for each 
unit. The estimates we will obtain on a will depend on the choice of norm. 
We choose: 

= max(|w|, \w\) , 

and consequently the matrix norm is chosen as: 
||M||= supY,\Mij\ . 



With this choice of norm we can write explicitly (for i = or 1): 

II inf , I,-"""'":'', , „ , (5.6) 

neSi neSi mecK[\ny,\ + \gw\,\ny\ + \gy\) 



with A{vn,Wn) — detA„. Now, this expression must be bounded in order to 
achieve an inequality of the form (3.8). The idea is to use the mean value 
theorem for the numerator and for the denominator of eq. (5.6). 

First, from eq. (5.1) it is seen that A(v^'^\ w^'^^) — A'^'^^ is positive. Hence, taking 
into account the condition on eq. (3.7), A{vn,Wn) must remain positive over 
the range of a value considered. Therefore we can drop the absolute value of A 
in expression (5.6). Now, A can be expanded as a Taylor series in its variables 
such that 

A{vn, wn, ) = A^"^ + {vn - v^^) + A« {w^ - «;«) + i?« , (5.7) 



where R^^^ is a remainder of order 2. We shall assume that i?^') is positive, 
which must be checked for given functions h and g. In particular, it is true for 
the FitzHugh-Nagumo model [eq. (5.3)]. (If R^^^ were not positive, we would 
have to estimate it in terms of (.) Then, using the variable ( introduced by 
eq. (3.1) we can bound A{vn,Wn) by 

A{vn,Wn)>A^'^-B^'^C, (5.8) 



where S« = + 

Second, the denominator of eq. (5.6) 
C is substituted such that we obtain 



is also expanded as a power series, and 
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\K\ + \g.\<C?+Dfc, 
\K\ + \g^\<C^'^ ^Dfc, 

where = \f^)\ + \g% = \h^\ + \g^\, and = |/^W| + + |/,W | + 
I g^l I , Df = \hfJ + \g'il\ + \h^^i\ + \ g^ I . These inequahties hold if the absolute 
value of the partial second derivatives |/i„t,(xW+0Ax)|, |/i„^(a;*^*)+^Aa;)|, etc. . ., 
with e e [0, 1] and Ax decreasing functions of 9. Again, this 

property is verified in the particular case of the FitzHugh-Nagumo model. We 
introduce the notation 

CW + D^^C = max(cP + D^C, + DfC) ■ (5-9) 



Finally, the substitution of eq. (5.8)-(5.9) into eq. (5.5)-(5.6) leads to 

II ||-i> min 7^,,^ . (5.10) 

We could insert this estimate into (3.7), but it would lead to a more compli- 
cated form of equation than (3.11). So instead we will bound (5.10) by 

||A-^||-i>mm(a«-a«C), (5.11) 

1=0,1 



where af = A^^/C^^,af = (>1«L>W +S«C(^))/C«2, and hi = {PC^^ - 
A^^^ Q) / C^^^^ . The reader can check this by multiplying both expressions by 

At this stage, we have obtained the inequality (3.8). In order to achieve the 
desired inequality (3.11) we have yet to determine the factors bo,Co and Ci 
defined by eq. (3.9)-(3.10). In the case of discrete diffusion, these factors are 
simply given by the following inequalities. 

2d 

\K{X) I = sup I ^ D{xi^ (a) - Xnia))\ 

neS 

<2ci||D|| sup Ixmict) — Xn{cy)\ 

m,n€S 

<2d\\D\\{C + 2C), (5.12) 
where ^ is equal to \x^^^ — x^^^\. 

Finally, for the discrete diffusion operator (2.3) it can readily be verified that 
the norm of DK is given by 

\\DK\\=4d\\D\\ . (5.13) 
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In conclusion, for a diffusive network of bistable elements described by two 
variables per unit we obtain the differential inequality 



da 



< 



2d \\D\\ (e + 2C) 
~ fliC — 4:d \\D II |q;| 



(5.14) 



which is of the form (3.11). Thus we can use (3.12) and obtain 



The FitzHugh-Nagumo system 

Let us deduce from (5.15) a lower bound of the critical coupling for which 
steady states persist in the network for the FitzHugh-Nagumo model defined 
by eq. (5.2)-(5.3). It is convenient to divide the second equation of the sys- 
tem (5.2) by 6, and this leads to better estimates. In this case the various 
factors of the inequality (5.10) are easily calculated and we can write 

II A-^ ||-i> min ( -'^P^ - '^\P^\C 

^=0'i Vmax(5 + e^y, jpl'V + \Pvv\CS + e) 



(5.16) 



Let us treat the case where the first term in the denominator is the largest 
one (which will be the case for all ( of interest if 7 is large enough) . Following 
eq. (3.8) we bound (5.16) by 

||A"^||"^> ao-aiC, (5.17) 



where oq — ai( is defined by mini=o,i(— pi'^ — \Pvv\0^'y / + ^7)- 

Now, II D 11= 1. The remaining constants to insert in (5.15) are ao, ai and ^. 
These depend on the form of the function p{v). As a concrete application, we 
consider the case where p{v) is a cubic polynomial defined by eq. (4.10). Then 
the equilibrium states of the bistable unit are at (0, 0) and (^, ^/^y) with 



Furthermore, in the present case we have p^^^ = —ri,p\,^' = r]— l-P^^j = 2(?7 + l) 
and p2^ = —2(2 — r/). So, if 77 < |, it can be shown that eq. (5.17) can written 
with 

_ e7?7 _ 2(?7 -I- l)e7 

"0 ~ 'x~, ' ^1 ~ — x~, ' 

+ e'-f + e'-f 
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Introducing the parameters a^, ai in eq. (5.15) we obtain that for the FitzHugh- 
Nagumo model continuation from the uncoupled limit is possible for ct < ckq 
with 



+ 0{ 




ao = 



8d{S + ej){r) + l)C 



We note that in the case 7 — cxd one recovers, at lowest order, the value 
ao ~ r]'^/8 which we obtained in eq. (4.11) for the one-variable cubic Nagumo 
model. 



6 Optical bistability 

In this section we will apply the scheme proposed in section 3 to an exam- 
ple taken from nonlinear optics. As mentioned in the introduction, one can 

conceive optical devices which, when submitted to an external driving field, 
may exhibit two different steady states. The first stable stationary state is 
characterized by a low transmission rate whereas the second stationary state 
has a high transmission rate. Thus, endowed with its two states, each device 
can cope with a unit of information, i.e., a bit. Therefore, a simple idea is 
to consider a network of such optical devices in order to create an optical 
memory. However, one of the problems which could arise in this context is 
that when the optical devices are brought close to each other, they become 
coupled, at least by means of their overlapping evanescent fields. The coupling 
might prevent the network from having stable stationary patterns. In this con- 
text, we propose in the present section to show how a bound on the coupling 
strength can be estimated, which guarantees the persistence of the stationary 
patterns which would exist in the absence of coupling. We note that this issue 
has been addressed by Firth in ref.[9], using a method closely related to the 
one of Keener [8] . 

As a simple model, we consider a network of identical optical ring cavities 
containing a nonlinear absorptive medium, and coupled in a d-dimensional 
lattice by means of their field amplitude. We assume that the dynamics of the 
network is governed by the following equations: 

Vn = -Vn + y - ^CWn + OL (f „+i - 2Vn + Vn-\) , 



= —{-Zn + 1 - VnWn) , 
K 

with n & Z. For simplicity the equations have been written with d — 1. 




(6.1) 
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In the absence of coupling (a = 0), the cq. (6.1) represent a simple three- 
variable model for studying absorptive optical bistability [14]. This model 
can be derived from the Maxwell-Bloch equations in a semi-classical mean- 
field approximation applied to a set of broadened two-level atoms [14]. The 
variables Vn,Wn, and reduced variables corresponding respectively to 

the real field amplitude produced in the cavity, the atomic polarization and 
the population difference between the lower and upper atomic levels. The 
model is also characterized by five parameters amongst which k, 7^ and 7|| are 
various decay rates. The amplitude of the field injected into the ring-cavity 
is represented by y and C corresponds to a combination of several physical 
parameters and is called the bistability parameter. For more details, see e.g. 
ref [14]. 

Although the bistable units are described by three variables, when these are 
coupled only through the variable v, the equilibrium states of the network can 
be studied by solving a one- variable system of the form (4.1) with 

fM = -v^ + y-2C-^. 

l + < 

Indeed, in the case where only v is coupled, the stationary states of and Zn 
can be expressed as function of Vn, namely w = f (1+w^)^^ and z = {l + v'^)~^. 
Therefore the scheme of section 4 apply. First we analyse the roots of f{v). In 
order to deal with analytical expressions for the stationary states we choose 
the value of parameter y as 

y^C+1. 

Then the stable equilibrium states of / are determined as 

= 1 , (6.2) 
= ^ (C + \/C2 - 4C - 4) ~ C - 1 . (6.3) 

Second we analyse the critical points of function /, i.e., the roots of f'{v). 
These are 

c± = ^C-l±^C{C-A). 

Now, as suggested in ref. [14], it is interesting to consider the limit of "fully 
developed hysteresis cycle" , that is the limit of large values of C. In this case 
the nearest critical points to {v^^\v^^^ is c_ because |c_ — v^^^ = 0{C^^) 
whereas |c+ — v^^^\ — 0{C). However, as |/'(w*-°'')| > [/'('w''^-*)!, we are in 
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a situation where there is an interval [0, Cm] about which \f'{v^^^ + Cm) I > 
_ Cm)|. This case is analogous to the one described by eq. (4.8). Hence 
following eq. (4.9) -or a slight generalisation when (v^^\v^^^) is not normalized 
to (0, 1)- continuation of the stationary states from the uncoupled limit is 
possible for a < ao where 

«0 = ^^^(,/_^(o)) (-/(^-) + /(^^°^ + Cm) + f{v^'^ - Cm)) . (6.4) 

The analytical expression of /(c_) is a complicated function of C we can 
expand as a power series in and we get 

f{c-)^~ + 0{C-'). 

On the other hand, the terms in (6.4) including Cm can be shown to be of 
order C~^. Hence, at the lowest order in C^^, eq. (6.4) gives 

«o = ^ + 0{C-') . 

This value of the coupling strength between the optical bistable cavities pro- 
vides a lower bound up to which the stationary patterns (2.4) continue to 
exist. Therefore this value of ccq gives a bound for which the optical network 
has a family of stationary states representing arbitrary indexed sets of binary 
digits. This is a prerequisite for designing optical memories. 

Finally note that we could consider diffusion of more than one variable in 
system (6.1). In this multivariable problem we could apply method as in sec- 
tion 3-5. 



7 Stability of the steady states 

In this section wc prove the perhaps obvious result that the continuation of a 
stable steady state of the uncoupled network is stable for small a. A word of 
caution is in order, however: it is not necessarily true that the stability-type 
is preserved for all a to which a steady state can be continued. It is perfectly 
possible that the equilibrium might undergo a Hopf bifurcation at some a : 
the equilibrium persists but changes stability-type. 

We treat spectral stability only. For a steady state X{a), we write M(a) = 
DG{X{a),a). The idea is that since M{a) varies continuously with a, its 
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spectrum can not move too fast. To make it precise, we introduce the resolvent 
of M: 

Rx{M) = (XI - M)-i , 

defined for A e C \ spec M. Then, using the formula A"^ — — A~^{B — 
A)B''^, one proves that A ^ specM' for all M' such that 

||M'-M||<||i?A(M)||-^ . (7.1) 

Now, let 

r= sup ||i?A(M(0)) II . 

ReA>0 

This is easy to calculate, as R\{M{0)) is block-diagonal, and by our hypotheses 
on spec DF(x*^*)) and boundedness of M, r < oo. Then it follows from (7.1) 
that M{a) has no spectrum in the right half plane (nor on the imaginary axis) 
as long as 

||M(q;) -M(0)||< J . 

Note that we can also continue unstable steady states of the uncoupled network 
for small a. They remain unstable for small enough a, because by the above, 
the resolvent remain bounded on the semicircle F with diameter (— i-B, +iB) 
on the imaginary axis, lying in the closed right half plane, for small a, where 
B is any number such that 

||M(0)||< S. 

Then there must continue to be spectrum inside the semicircle F because the 
spectral projection operator Pr — Ir R\{M{a))dX depends continuously on 
a and is non-zero for a = 0. 



8 Conclusion 

The aim of this paper was to study the stationary states of networks of bistable 
units with weak coupling. We proved that, for a large class of network inter- 
action, not necessarily linear, and for small values of the coupling strength, 
the network presents a high steady state multistability inherited from the 
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bistable behaviour of each unit. In some cases, it is possible to give a boTind 
on the coupling strength up to which the stationary solutions continue to 
exist. We have applied this analysis to concrete examples, such as the two- 
variable FitzHugh-Nagumo system and a simple model describing an instance 
of optical bistabihty. 

It is plausible that the existence of so many stable stationary states would 
block wavefronts from propagating in the network. This phenomenon is called 
'propagation failure', and has been shown to be an important discreteness 
effect in diffusion-coupled bistable systems [8]- [12]. We wish to emphasize that 
our present results, however, do not say anything about the existence or not of 
propagating waves or fronts. Nevertheless, we believe that our work could be 
extended to prove structural stability of the uncoupled case, i.e. topological 
equivalence of the weakly coupled case to the uncoupled case. This would 
preclude existence of propagating waves of interest, as the uncoupled case 
has no propagating wave solutions, except those involving the unstable steady 
state between the two stable steady states. These are probably not of any 
interest as they require preparation of the initial condition in a multiply- 
unstable situation. 

The present study illustrates a simple and efficient method for investigating 
weakly coupled systems. The method consists in analysing, first of all, the 
solutions of the system in the limit of vanishing coupling in the network. 
This limit is interesting because it may provide an infinite number of explicit 
solutions each of which having possibly a unique continuation for non-zero 
value of the coupling. Moreover, it is possible, in principle, to estimate the 
range of coupling strength over which the continuation from the uncoupled 
limit exists. 

Note that, as in [7,13], we could also deduce a locahsation length for the steady 
states for small a though we have not done so here. This means that the effect 
of changing the state of one unit decays exponentially from that site. For small 
a, and local coupling, the decay rate is 0{a). 

The procedure proposed in this study could probably be applied to more com- 
plex situations than stationary states. For example, the units of the network 
could be replaced by bistable oscillators which have a stable stationary state 
in addition to a stable limit-cycle. In that case it is interesting to determine 
whether there exist stable localized oscillations, as a continuation from the 
uncoupled limit. In fact this question has already been addressed and an- 
swered affirmatively in the case of Hamiltonian or time-reversible oscillatory 
networks [13]. For such systems, it has been recently proved that there exist 
'discrete breathers', i.e., localized oscillatory solutions in the network, when 
the coupling strength between the oscillators is sufficiently small. Notice that 
in that case, the operator G defined by eq. (2.5) corresponds to an operator de- 
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fined on 'spaces of T-periodic loops', in order to cope with the time-dependent 
behaviour of the oscillators. For dissipativc networks, the sort of results we 
could expect would be a network analogue of the theory of persistence of 
normally hyperbolic invariant manifolds. 
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